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We calculate the phase, the temperature and the junction length dependence of the supercurrent 
for ballistic graphene Josephson-junctions. For low temperatures we find non-sinusoidal dependence 
of the supercurrent on the superconductor phase difference for both short and long junctions. The 
skewness, which characterizes the deviation of the current-phase relation from a simple sinusoidal 
one, shows a linear dependence on the critical current for small currents. We discuss the similarities 
and differences with respect to the classical theory of Josephson junctions, where the weak link is 
formed by a diffusive or ballistic metal. The relation to other recent theoretical results on graphene 
Josephson junctions is pointed out and the possible experimental relevance of our work is considered 
as well. 
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I. INTRODUCTION 

The peculiar electronic properties of graphene first ob- 
served experimentally by Geim et alji and Zhang et al.^ 
can accurately be described by massless two dimensional 
Dirac fermion excitations (for reviews on the physics of 
graphene see, e.g, Refs. 

Owing to the proximity effect a superconductor can 
induce non-zero pair-potential in the graphene as well^. 
Such graphene-superconductor hybrid structures and in 
particular the Andreev-reflection taking place at the 
graphene-superconductor interface was first studied the- 
oretically by Beenakker* (for a review on Andreev re- 
flection in graphene see Ref. 0). Soon after Beenakker's 
pioneering work, supercurrent between two supercon- 
ducting electrodes on top of a graphene monolayer has 
been observed exp erimentally by Heersche et al lid , and 
later in Refs. IllMlSl In particular, the experimental 
results of Ref. |T2 attest to the ballistic propagation of 
quasiparticles in graphene-superconductor hybrid struc- 
tures, whereas the experiment of Du et al}'^ gave evi- 
dence that it was possible to fabricate transparent SG 
interfaces. These experiments have also sparked consid- 
erable theoretical interest in superconductor-graphene- 
superconductor (SGS) heterostructures. The short junc- 
tion limit, where the coherence length ^ = hvp /Aq (here 
Di? is the graphene Fermi velocity and Aq is the supercon- 
ducting gap) is smaller than the length L of the junction, 
was first studied by Titov and Beenakker— assuming 
ballistic graphene. In the opposite, long junction limit 
the density of states of the Andreev levels was calcu- 
lated first by Titov, Ossipov and Beenakker—. Subse- 
quently, numerous other theoretical works investigated 
the Josephson current in SGS structure o"'^^'^^' '— . The 
tunneling effect in SG structures has been studied in 
several works^^"— as well. Other works in the field of 
graphene-superconductor heterostructures include stud- 
ies on crossed Andreev reflection in a graphene bipolar 
transisto r^^i'^^ , on s- and d-wave SG junctions^iH and on 
ferromagnetic SG structures^i^. Very recently, using a 



phase-sensitive SQUID interferometry technique Chialvo 
et al. has studied experimentally the current-phase rela- 
tion (C'I'R) of graphene Josephson junctionsiS. 

In this work we calculate the Josephson current in 
SGS structure as a function of the superconductor phase 
difference, the temperature and the length of the junc- 
tion. In our theoretical treatment we adapted the method 
used by Brouwer and Beenakker for metallic chaotic 
Josephson junctions^. The approach allows to obtain 
results for finite temperature and is valid for junctions 
of arbitrary length. We note that this method has al- 
ready been applied for calculating the persistent current 
through a n-p junction in graphene^^. Wherever possi- 
ble, we compare our results to previous ones derived for 
superconductor-normal conductor-superconductor (SNS) 
junctions, where the normal conductor is a ballistic 
metal. 

The rest of the paper is organized in the following way: 
in the next section we introduce the theoretical approach 
that we use to obtain the Josephson current. In Sec- 
tion mil we present and discuss the results of numerical 
calculations. Finally, in Section llVI we give a brief sum- 
mary. 



II. THEORETICAL CONCEPT 

We consider a Josephson junction in the x-y plane. 
The normal graphene region (G) at < L/2 separates 
the two superconducting regions formed by covering the 
graphene layer by two superconducting electrodes (S) in 
the regions x < —L/2 and x > L/2 (for the geome- 
try see Ref. [TtI ). The width of the Josephson junction 
along the y axis is W. Owing to the valley degeneracy 
of the Hamiltonian, the full Dirac-Boguliubov-de Gennes 
(DBdG) equations for graphene-superconductor systems 
decouple to two four by four, reduced Hamiltonians that 
are related to each other by a unitary transformation 
(see, e.g., Ref. HI). We now take the one corresponding 
to the valley K. Then the quasi particle excitations in 
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the SGS systems are described by the reduced DBdG 
equations: 



Ho-fi A(x,y) 



(1) 



where iJo = —ihvF{crxdx + crydy) + U{x,y) ao is the Dirac 
Hamiltonian. Here ax and ay are Pauli matrices, ctq is the 
unit matrix and, /x is the chemical potential and £ > is 
the excitation energy. The superconductor electrodes are 
doped by the potential U{x,y) = UqQ{\x\ — L/2) (here 
Uq < and constant, and Q{x) is the Heaviside func- 
tion). The wave function = (^e,5';j) is comprised 
of electron and hole wave functions which have 
opposite spin and valley indices. For the pair potential 
A(a;, y) we assume a simple model: its magnitude Ag 
is constant in the S regions, changes step-function-like 
at the SG interfaces (so called "rigid boundary condi- 
tion", see Ref. 37) and is zero in the normal conducting 
region. Similarly, we assume that the its phase is piece- 
wise constant in the S regions. Hence, the pair poten- 
tial is given by A{x,y) — Aq e~"^/^ for x < —L/2 and 
A{x,y) = Aoe*'^/^ for x > L/2. Band bending or other 
effects of the sup erconducting electrodes are neglected 
(see e.g. Ref. |33 for the discussions of some of these ef- 
fects in the case of normal conducting metal electrodes). 

The Josephson current at finite temperature is given 
by26 



I=-2k^T'4i- 



deg{e) In 



2 cosh 



2kBT 



(2) 

where (j) is the phase difference across the junction, g{e) 
is the density of states of the Andreev levels. The factor 
of 4 accounts for the spin and valley degeneracy. As 
one can see from Eq. a necessary step to calculate 
the Josephson current is to obtain the density of states 
of Andreev bound states in the SGS junction. To this 
end one can in principle proceed in the following way: 
one can write down a trial wave function in all three 
regions of the SGS structure. The boundary conditions 
at the two graphene-superconductor boundaries of the 
SGS junction then result in a secular equation J-{e) — 
whose solutions ei give the energies of the Andreev bound 
states. For finite Uq we obtained a 8 x 8 determinant 
for the secular equation (since this determinant is quite 
lengthy here we do not present its detailed form). Once 
the energy levels Si of the SGS junction are known, the 
density of state is given by g{e) = J2i S{e — Si). 

However, the above outlined method is numerically 
quite cumbersome since one has to search for the ze- 
ros of the secular equation J^{e) = 0. To overcome this 
problem we now follow the method used by Brouwer and 
Beenakker in Ref. H^. They rewrote the expression for 
the Josephson current given in Eq. ([2]) in a more conve- 
nient form. Here we only summarize the main steps of 
the derivation. The secular equation can be written as 
J-{e) = ToIli{{e — Si)) = 0, where Tq is a function of e 



but does not have zeros in the complex plane. Thus it is 
easy to see that the density of states can be expressed as 



g{e) 



1 d 

TT de 



Im In J'(e + iO^ 
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where 0^ is a positive infinitesimal. Using the analytic 
properties of J- in the upper half of the complex e-plane 
and the fact that under the change e ^ —e the func- 
tion T goes over into its complex conjugate (physically, 
this follows from the electron- hole symmetry), the e- 
integration can be extended from — oo to oo. Finally, 
after performing a partial integration the Josephson cur- 
rent in Eq. ^ can be rewritten as 



^ _ 2e d 
inh d(j) 



oo+iO+ 

de tanh 

-oo+iO+ V^'SbJ 



In^(e). (4) 



Now closing the integration contour in the upper half of 
the complex e-plane and applying the residue theorem 
the Josephson current becomes 



4 d 

n=0 



h 



(5) 



where ia;„ = i(2n + l)7rfcB2^ are Matsubara frequencies. 
Note that in our model lnJ^(£) has no singularities for 
Im £ > 0, thus the poles of the integrand in Eq. ^ 
come only from the hyperbolic tangent function. The 
main advantage of this result is that one does not need 
to obtain explicitly the solutions of the secular equation 
J-{e) = 0. Moreover, this method immediately gives the 
finite temperature dependence of the Josephson current. 
In the numerical calculations it turns out that the sum 
in Eq. ([5]) is rapidly convergent and usually one need to 
include only a finite number of terms. A similar result 
has been found for the persistent current through a n-p 
junction in grapheneSS.. 

We now consider the experimentally relevant case of 
highly doped superconductor electrodes, ie, the limit 
Uq — > — oo. By matching the wave functions at the 
graphene-superconductor boundaries of the SGS struc- 
ture we found the same secular equation J-"(£, q-m) = as 
that obtained by Titov and Beenakker using the transfer 
matrix method (see Eq. (14) in Ref. [13) ■ We used the 'in- 
finite mass' boundary conditions'^^ at y = and y = W 
for which q„i — {m + l/2)Tr/W, where m — 0, 1, 2, . . . 
(for W ^ L the choice of the boundary conditions is ir- 
relevant). For a given m the solutions of the quantization 
condition J-{s, Qm) = give the Andreev energy levels £„ 
for £,„ < |A|. The secular equation T{e, g™) = is valid 
both in short and long junction limiti^. One can show 
that J-{—e + iO'^, q) = T*{e + iO~^, q) which is a necessary 
condition"^^ for writing the Josephson current in the form 
of Eq. Q. 

The current contribution from each propagating mode 
with transverse momentum q^ can be calculated sepa- 
rately and the total current is the sum over these contri- 
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butions: 

^ M oo ^ 

I = ~y ^fceT X! X! ^ \nT{iuJn,qyn), (6) 

m=l n=Q ^ 

where M is the number of propagating modes and the 
function J-{e,q) determines the energy levels for a given 
transverse momentum q. This equation is our starting 
point for calculating the supercurrent through a graphene 
based Josephson junction. 

Further analytical progress can be made in the short 
junction limit (L <C because the Andreev levels can 
in that case be obtained in a closed form (see Eq. (16) 
in Ref. Il7t ). Similarly, the summation over the Mat- 
subara frequencies in Eq. ([B]) can be performed ana- 
lytically using the identity T^.T ^ 1/ [(2fc + 1)^ + x^] = 
(7rtanh(7ra::/2)/(4a;) (see Ref. |40|) and we find 




Here e™ and Tm can be found in Ref. [itI while the tem- 
perature dependence of the superconductor gap |A| — 
Aq(T) for s-wave superconductors is given by 

^ ^ n— 1 ^ ^ 

where Kq{x) is the zero order modified Bessel function, 
Ao(0) = (eT/7r)fcBrc = O.SGTfcB?; and 7 is the Euler's 
constant^. For zero temperature from ^ one can ar- 
rive at the same expression for the Josephson current as 
that obtained by Titov and Beenakker (Ref. IT7| ). For 
finite temperatures the summation over the transversal 
modes m in Eq. (jB]) cannot be evaluated analytically but 
numerically can easily be treated. 



|i/Ao(0) = Q |a/Aq(0) = 20 
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FIG. 1. (color online) The supercurrent (in units of Jo = 
eAo(0)/fi) as a function of the phase difference </!>. The pa- 
rameters are as follows: In (a) and (b) [short junction limit] 
^/L = 20 and T/Tc = 0, 0.53, 0.71 (black □, red o and blue A, 
respectively). In (c) and (d) [here L ~ ^] we used ^/L = 0.91 
and T/Tc = 0,0.18,0.35 (black □, red o and blue A, respec- 
tively). In (e) and (f) we consider long junctions [^/L = 0.05]. 
In (e) T/Tc = 0,0.035,0.053 (black □, red o and blue A, re- 
spectively). In (f) T/Tc = 0,0.018,0.035 (black □, red o and 
blue A, respectively). The chemical potential is /x = in (a), 
(c) and (e), whereas it is /i/Ao = 20 in (b), (d) and (f). The 
width of the junction is £,/W = 0.05 in all cases. 



III. NUMERICAL RESULTS 

We now present the results of numerical calculations 
for the Josephson current using the most general for- 
mula given by Eq. (jG]). Figure [1] shows the supercur- 
rent as a function of phase (f) for a number of interesting 
case: for short (L <C ^), intermediate (i ~ ^) and long 
(L S> C) junctions, assuming ^ — and for finite /i/Ap as 
well. One can make the following general observations: 
a) the maximum current increases by increasing the dop- 
ing (p/Aq value) and by decreasing the temperature or 
the junction length; b) at higher temperatures the cur- 
rent shows a simple sinusoidal dependence on the phase 
in all of the cases while at low temperatures the position 
of the maxima of the currents are shifted to the right 
resulting in a skewness of the curves. Following Ref. [l6| 
we define the skewness by = 2(/),nax/7r — 1, where ^max 
is the position of the maxima of the supercurrent. Both 
the tendency to simple harmonic dependence for T — >■ Tc 
and a positive skewness (i.e. i^max > t^/^) are in line 



with previous calculations on Josephson current in weak 
links which comprise a normal conducting metal or a tun- 
nel barrier and assume that the pair potential changes 
abruptly at the normal-superconductor interface (see e.g. 
Ref. I33 for a review) . 

We find especially interesting the results shown in 
Figs. [TJe) and (f). In the long junction limit for /i = 
we find that the skewness is very small even at T/Tc — 
[the curve denoted by black squares in Fig. [T{e)] thus 
the C$R resembles a harmonic dependence. In contrast, 
still in the long junction limit but for /i/Ao = 20 and 
T/Tc = 00 [black squares in Fig. [Ijf)] we see that the 
current depends almost linearly on the phase and the 
curve resembles a saw-tooth. (The transition to a saw- 
tooth-like dependence can already be seen in Fig. [IJd) 
where L/^ = 1.1, c.f. Fig. [IJb) showing the short junc- 
tion limit.) It is interesting to note that the theoretical 
result for clean, long SNS junctions at low temperature is 
a saw-toothed C^R^i^. Our numerics suggests that for 
SGS junctions in the same limit the saw-tooth is some- 
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what rounded-ofF and the slope of the curve is finite when 
(/) TT. Thus, the C$R in long, clean SGS junction seems 
to be closely resembling of, but not identical to the cor- 
responding case in SNS junctions. 
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The temperature dependence of the pair potential was 
taken into account using Eq. The results are shown 
in Fig. [3] At this point it is interesting to make a quick 
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FIG. 2. (color online) The skewness S as a function of the crit- 
ical current Ic for different coherence lengths ^. The parame- 
ters are ^/L = 0.35, ^/W = 0.0077 (black dots), ^/L = 1.05, 
^/W = 0.0231 (red squares) and ^/L = 1.75, ^/W = 0.0385 
(blue triangles). The lines are guides to the eye. The ratio 
of the chemical potential and the superconducting gap was 
fi/Ao = 10. The dotted line shows that for small critical 
currents the skewness depends linearly on Ic- 

We calculated the skewness 5 as a function the crit- 
ical supercurrent Ic (the value of the current at 0max) 
and plotted the results for three different ^ values, while 
keeping the junction length L and width W constant. 
The used ^/L values go from ^/L — 0.35 (long junction 
limit) to ^/L ~ 1.75 (short junction limit). There are 
two important things to notice in Fig. [2] for small criti- 
cal currents: a) for a given junction length L, as /c — > 
(for higher temperatures) the skewness S also goes to 
zero, i.e. the C<i>R is approaching a simple sinusoidal 
form; b) S depends linearly on 1^ for small critical cur- 
rents, while at larger Ic the dependence clearly deviates 
from a simple linear relation. The skewness has recently 
been measured in Refd and the case of ^/L = 0.35, 
^/VF = 0.0077, shown by black dots in Fig. [21 in princi- 
ple corresponds to that of sample B in this experiment 
(with estimated coherence length of ^ w 100 nni^). Our 
calculations give a smaller slope than the measurements 
in Ref. [H. According to our numerics, the larger slopes 
observed in this experiment would be attainable in the 
short junction limit. Note however, that the exact slope 
would also depend on the value of the chemical potential 
which was not known, and importantly, the samples in 
the experiment of Ref. JJj; are likely to have been in the 
quasi-diffusive limit, therefore we cannot expect quanti- 
tative agreement with our ballistic theory. 

We have also calculated the temperature dependence 
of the critical current for short, intermediate and long 
junctions, taking two values of the chemical potential /i. 
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FIG. 3. (color online) The critical current Ic as a function of 
TjTc. The parameters are ^/L = 20 in (a), ^jh = 0.91 in (b) 
and ^jL = 0.05 in (c). In (b) and (c) we used logarithmic 
scale. Red o denote the results for p = and black □ for 
/i/Ao(0) = 20. The width of the sample was fixed : ijW = 
0.05. 



sidestep and note that Titov and Beenakker (Ref. Ill 
showed that for short junctions, at zero temperature 
and at the Dirac point the C$R for ballistic graphene 
is formally identical to the classical result of Kulik and 
Omel'yanchuk^, which, however, assumes diffusive nor- 
mal metal as a weak link. Looking at the Ic — T curves in 
the short junction limit [ Fig.[21[a)] one can see that they 
are also qualitatively similar to the corresponding result 
of Kulik and Omel'yanchuk^ [c.f. Fig. 7 in Ref. Is^. 
The close resemblance of certain properties the two types 
of Josephson junctions therefore seems to extend to the 
temperature dependence of the critical current as well. In 
the opposite, long junction regime, for ^ — Q [shown by 
red circles in Fig.[3l^c)] one can observe a short plateau in 
the current for small temperatures followed by an expo- 
nential decay. Interestingly, a qualitatively very similar 
result has been obtained by Gonzalez and Perfetto in 
Ref. [23l . assuming tunnel coupling between the supercon- 
ductor and the graphene and using a different formalism 
than ours. In the case of finite doping [black squares in 
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Fig. [3jc)] an exponential decay of Ic can be seen basi- 
cally in the whole temperature range. The 1^ in clean 
SNS junction exhibits the same qualitative dependence 
on T (Ref.il). 

Experimentally, the Ic ~ T relation was measured by 
Du et al. (Ref. and by Ojeda-Aristizabal et al. 
(Ref. 15). Again, quantitative comparison with our re- 
sults is not possible because the graphene samples in the 
experiments were in the diffusive limit, but the observed 
dependence of I,, on the temperature was qualitatively 
similar to the results shown in Fig. ^a). 

Finally, we studied the length dependence of the crit- 
ical current and the results are shown in Fig. S) At the 



is remarkably close to the results of the self-consistent 
tight-binding calculations of Black-Schaffer and Doniach 
(Ref. [13) who found b = —1.3. Considering now the 
case of doped graphene weak link, we found that for 
L/^ > 1 the L dependence of Ic can be well fitted by 
Ic/Iq = laC^^^^^ with 6 ss 1. Exponentially small criti- 
cal current is also typical for clean SNS junctions if L is 
larger than the thermal coherence length (for details see 
e.g. Ref. d?!)- We note that in the case of L/^ < 1 and 
T/Tc = one can see oscillations in the Ic vs L curve 
[shown in the inset of Figlljb)] whose study is left to a 
future work. 




FIG. 4. (color online) The critical current Ic as a function of 
the junction length L/^ in logarithmic scale. In (a) we used 
/i = 0. Black □ and red o denote the results of T/Tc = 0.0 
and T/Tc = 0.06 calculations, respectively. The symbols in 
the inset of (a) show Ic for T/Tc = 0.0, L/^ < 1 in double 
logarithmic plot, along with the fitted linear function (solid 
line, see text). In (b) the chemical potential is /i/Ao(0) = 10. 
Black □ and red o denote the results of T/Tc = 0.0 and 
T/Tc = 0.18 calculations, respectively. The inset of (b) shows 
the T/Tc = 0.0 calculations in linear scale for L/^ < 2. The 
width of the junction was £,/W = 0.05 in all cases. 

Dirac point {fi = 0) one can observe an exponential de- 
cay of Ic for L/^ » 1 [main panel of FigUJa)]. We could 
not see the ~ 1/L^ dependence predicted in Ref. [1^ For 
■^/? ^ 1 and T/Tc = 0, however, we do find a power- 
law dependence Ic ^ T^ [see the inset of FigUJ^a)] and 
fitting the numerical results we obtained h = —1.4. This 



IV. CONCLUSIONS 

In this work we calculated the Josephson current in 
ballistic SGS structures. The most important assump- 
tion we made is that one can use rigid boundary con- 
ditions, i.e. that the change of the pair potential is 
step-function-like at the SG interfaces and that it does 
not depend on the supercurrent. We developed a gen- 
eral and numerically efficient approach to obtain the 
current-phase relation for arbitrary length of the junc- 
tion as well as for finite doping and temperature. At 
low temperatures we have found that the current-phase 
relation differs from a simple harmonic dependence. In 
the case of short junctions and small critical currents 
the deviation of the current-phase curve from the sinu- 
soidal form, the skewness, shows a linear dependence on 
the critical current, similarly to the observation of a re- 
cent experimentiSj though the slope of the curve did not 
match the experimental one. This is likely to be due 
to the fact that in the experiment the graphene sample 
was quasi-diffusive. In the long junction limit our results 
show that in clean SGS junctions the the current-phase 
relation transforms from the sinusoidal form at T < Tc to 
a curve resembling saw-tooth at T ^ Tc. In contrast to 
clean SNS junctions however, our numerics suggests that 
the dependence is not exactly saw-toothed. We have also 
calculated the temperature and junction length depen- 
dence of the critical current. We have found similarities 
to both classical results for SNS junctions and recent ones 
obtained for graphene but using a different formalism^^. 
In respect of these numerical calculations further theo- 
retical progress is needed to unravel the relation between 
the SGS and SNS results. 

Since the fabrication of both ballistic graphene samples 
and transparent SG interfaces have already been demon- 
strated experimentally^^'^14, we believe that our theoret- 
ical approach may be useful in the future to understand 
and analyze experimental data. In particular, the mea- 
surement of the C$R and the length and temperature 
dependence of Ic in short junctions should be feasible. 

Note added: During the peer- review process of the 
manuscript, a relevant preprint has appeared*^^ where 
the authors study the temperature dependence of the 
Josephson current using self-consistent tight-binding nu- 
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merical computations. Their work is thus complementary 
to our and help to understand the scope of certain ap- 
proximations, e.g. the rigid-boundary condition that we 
employed. 
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